#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include "matrix_utils.h"
#include <math.h>

#define LOWER 1e-9
// 求逆矩阵
double** matrixInverse(int *rows, int *cols) {
	char continueInput[100];
	double **matrixIn = NULL; 
	
	double matrixA(double **mar, int n);
	
	do{
		if (!check(rows, cols)) {
            continue;
        }
        if (*rows != *cols) {
            printf("此矩阵无法求逆，请重新输入！\n");
            continue;
        }
        
		double **newMatrix = enterMatrix(*rows, *cols);
    	if (newMatrix == NULL) {
        	return NULL;
    	}
    	double num = matrixA(newMatrix, *rows);
	    if (fabs(num) < LOWER) {
			printf("此矩阵不可逆，请重新点击选项输入！\n");
			freeMatrix(newMatrix, *rows);
			return NULL;
		}
	    // 创建增广矩阵 [A|I]
	    double **augmentedMatrix = (double **)malloc(*rows * sizeof(double *));
	    if (augmentedMatrix == NULL) {
			return NULL;
		} 
	    for (int i = 0; i < *rows; i++) {
	        augmentedMatrix[i] = (double *)malloc(2 * *cols * sizeof(double));
	        if (augmentedMatrix[i] == NULL) {
			for (int j = 0; j < i; j++) {
				free(augmentedMatrix[j]);
			}
			free(augmentedMatrix);
			freeMatrix(newMatrix, *rows);
			return NULL;
			}
	        for (int j = 0; j < *cols; j++) {
	            augmentedMatrix[i][j] = newMatrix[i][j];
	            augmentedMatrix[i][j + *cols] = (i == j) ? 1.0 : 0.0;
	    	}
	    }
	
	    // 高斯 - 约旦消元
	    for (int i = 0; i < *rows; i++) {
	        // 选主元
	        if (fabs(augmentedMatrix[i][i]) < LOWER) {
	            int pivotRow = -1;
	            for (int k = i + 1; k < *rows; k++) {
	                if (fabs(augmentedMatrix[k][i]) > LOWER) {
	                   	pivotRow = k;
	                    break;
	                }
	            }
	            if (pivotRow == -1) {
	                printf("矩阵不可逆。\n");
	                freeMatrix(augmentedMatrix, *rows);
	                continue;
	           	}
	            // 交换行
	            for (int j = 0; j < 2 * *cols; j++) {
	                double temp = augmentedMatrix[i][j];
	                augmentedMatrix[i][j] = augmentedMatrix[pivotRow][j];
	                augmentedMatrix[pivotRow][j] = temp;
	            }
	        }
		    // 将主元化为 1
			double pivot = augmentedMatrix[i][i];
		    for (int j = 0; j < 2 * *cols; j++) {
		        augmentedMatrix[i][j] /= pivot;
		    }
		    // 消去其他行的第 i 列元素
		    for (int k = 0; k < *rows; k++) {
		        if (k != i) {
		            double factor = augmentedMatrix[k][i];
		            for (int j = 0; j < 2 * *cols; j++) {
		                augmentedMatrix[k][j] -= factor * augmentedMatrix[i][j];
		            }
		        }
			}
		}
	
		// 提取逆矩阵
		matrixIn = (double **)malloc(*rows * sizeof(double *));
		if (matrixIn == NULL) {
			return NULL;
		} 
		for (int i = 0; i < *rows; i++) {
		    matrixIn[i] = (double *)malloc(*cols * sizeof(double));
		    if (matrixIn[i] == NULL) {
			for (int j = 0; j < i; j++) {
				free(matrixIn[j]);
			}
			free(matrixIn);
			free(augmentedMatrix);
			freeMatrix(newMatrix, *rows);
			return NULL;
			}
		    for (int j = 0; j < *cols; j++) {
		        matrixIn[i][j] = augmentedMatrix[i][j + *cols];
		    }
		}
		freeMatrix(augmentedMatrix, *rows);
		    
		printf("------原矩阵----------\n");
		printMatrix(newMatrix, *rows, *cols);
		printf("------运算后----------\n");
		printMatrix(matrixIn, *rows, *rows);
		
		freeMatrix(newMatrix, *rows);
		
		printf("是否继续输入矩阵（y(Y) 继续，其他值退出，谢谢合作！）\n");
		while (getchar() != '\n');
		fgets (continueInput, sizeof(continueInput), stdin);
        continueInput[strcspn(continueInput, "\n")] = 0;
    } while ((continueInput[0] == 'y' || continueInput[0] == 'Y') && continueInput[1] == '\0');
	
    return matrixIn;
}
// 求行列式的值 
double matrixA(double **mar, int n){
	double A = 0;
	double sign = 1;
	// 二阶方阵的情况
	if (n == 2) {
		return mar[0][0] * mar[1][1] - mar[0][1] * mar[1][0];
	} else {
		// 动态分配数组指针内存 
		double **newmar = (double**)malloc((n - 1) * sizeof(double*));
		if (newmar == NULL) {
			return 0;
		}
		// 动态分配行指针内存
		for (int i = 0; i < n - 1; i++) {
			newmar[i] = (double*)malloc((n - 1) * sizeof(double));
			if (newmar[i] == NULL) {
				for (int j = 0; j < i; j++) {
					free(newmar[j]);
				}
			free(newmar);
			return 0;
			} 
		}
		
		for (int j = 0; j < n; j++) {
			int ni = 0;
			for (int i = 1; i < n; i++) {
				int nj = 0;
				for (int k = 0; k < n; k++) {
					if (k == j) {
						continue;
					}
				newmar[ni][nj] = mar[i][k];
				nj++;
				}
			ni++;
			}
		A += sign * mar[0][j] * matrixA(newmar, n - 1);
		sign = -sign;	
		}
		// 释放子矩阵的内存
		for (int i = 0; i < n - 1; i++) {
			free(newmar[i]);
		}
		free(newmar);
	}
    while (getchar() != '\n');
	return A;
}

